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Abstract 


The mechanical response of proton exchange membranes in a fuel cell assembly is investigated under humidity cycles at a constant temperature 
(85°C). The behavior of the membrane under hydration—dehydration cycles is simulated by imposing a humidity gradient from the cathode to 
the anode. Linear elastic, plastic constitutive behavior with isotropic hardening and temperature and humidity dependent material properties are 
utilized in the simulations for the membrane. The evolution of the stresses and plastic deformation during the humidity cycles are determined 
using finite element analysis for two clamping methods and various levels of swelling anisotropy. The membrane response strongly depends on 
the swelling anisotropy where the stress amplitude decreases with increasing anisotropy. These results suggest that it may be possible to optimize 
a membrane with respect to swelling anisotropy to achieve better fatigue resistance, potentially enhancing the durability of fuel cell membranes. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Proton exchange membrane fuel cells (PEMFCs) are viable 
candidates for powering vehicles in a proposed future hydrogen 
economy. However, in order to replace the internal combustion 
engine, the durability and reliability of the fuel cell systems must 
be improved. For automotive applications the current target is 
5000 h (150,000 miles equivalent) operational life over a full 
range of vehicle operating temperatures (—40 to 40°C) [1]. 

During the operation of the cell, the proton exchange mem- 
brane (PEM) provides an ionic conductive path for protons from 
the anode to the cathode, while acting as an electronic insulator 
and a gas barrier to prevent mixing of oxygen and hydrogen. Any 
type of discontinuity in the membrane reduces the performance 
(i.e. output power, total voltage) and lifetime of the cell [2—6]. 
Although the electro-chemical and thermo-mechanical interac- 
tions among cell components (electrocatalysts, membranes, gas 
diffusion layers, and bipolar plates) affect the durability, it has 
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been found that the membrane itself is a major source of fail- 
ure, including mechanical damage and chemical degradation 
[2,3,5—8]. Thus, the membrane must be durable enough to with- 
stand mechanical stresses and chemical attacks so that the fuel 
cell can sustain its function under the severe internal operating 
conditions. 

Several forms of mechanical damage in the membrane elec- 
trode assembly (MEA) are commonly observed, including 
through-the-thickness tears or pinholes in the membrane, or 
delamination between the polymer membrane and the electrodes 
[4-10]. A number of studies related to the influence of the 
membranes’ mechanical properties on the durability and per- 
formance of the cell are available in literature [2,3,9,11-13]. 
It is commonly believed that the mechanical stresses, due to 
the hydration—dehydration cycles in the membrane precipitate 
these damage mechanisms [9,14—16]. Kolde et al. [13] investi- 
gated the material properties and fuel cell performance of several 
commercial fuel cell membranes of various thicknesses. The 
results indicated that membranes exhibiting good dimensional 
stability in the in-plane directions during hydration—dehydration 
cycles are desirable for improving the cell performance and 
product reliability. A key parameter is therefore the swelling 
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coefficient of the membrane [14,16,17]. Previous work sim- 
ulating the mechanical response of a PEM during fuel cell 
operation [14] showed (through numerical means) that compres- 
sive stresses develop upon hygro-thermal loading. In some cases, 
these stresses can exceed the yield strength, causing permanent 
deformation. This, in turn, results in tensile residual stresses after 
dehydration. These in-plane tensile residual stresses are believed 
to be a significant contributor to the mechanical failures observed 
in the membranes, since they may cause the propagation of the 
through-the-thickness cracks in the membrane. 

Even though crack initiation mechanisms in polymers are 
not completely understood, cyclic loading results in fatigue 
crack growth, which is one of the mechanisms governing the 
lifetime of these materials. The fatigue crack growth rate is gov- 
erned by the stress amplitude and stress level [18,19]. Therefore, 
investigating the stresses in the membrane during hygro-thermal 
cycling, which is the objective of this paper, is an important part 
of understanding the failure mechanisms in the membrane. 

A common method to experimentally simulate the fuel cell 
operation is with accelerated cycling through a range of inlet 
relative humidities (RH) at elevated temperatures (80-95 °C) 
(3,7,13,15,20,21]. Accelerated open-circuit voltage decay tests 
were conducted by Protsailo [20] using Nafion® membranes 
to identify the long-term effects of high temperature and low 
RH operation on the performance of the perfluorosulfonic 
acid (PFSA) membranes. They found that membrane lifetime 
increases with increased inlet humidity and decreases with an 
increase in operating temperature. In the case of reduced relative 
humidity, the membrane failure is accelerated by a degrada- 
tion of the mechanical properties of the PFSA [20]. In order to 
demonstrate purely mechanical failure modes, several investiga- 
tors conducted RH cyclic tests in the absence of electro-chemical 
reactions [7,15]. Mathias et al. [15] conducted an RH cycling test 
for an MEA made from a cast, 25 pm N afion®! membrane. The 
MEA was held at a constant 80 °C and cycled between maximum 
and a minimum RH values in the absence of chemical degrada- 
tion [15]. Recently, different test procedures were introduced 
to study the purely mechanical, and combined mechanical and 
chemical durability of fuel cell membranes for the development 
of new materials [7,21]. Crum and Liu [7] conducted accelerated 
and non-accelerated tests (representing realistic automotive duty 
cycles) for various GORE-SELECT®* membranes (produced 
by W.L. Gore & Associates Inc.). These tests were run at a con- 
stant temperature (80°C) and varying inlet and outlet RH at the 
anode and cathode. In order to separate the effects of chemical 
degradation and mechanical fatigue, they performed the tests in 
an inert environment (nitrogen) to suppress chemical activity. 
Similar trends in the durability of membranes were observed for 
the both accelerated and non-accelerated tests [7]. Thus, accel- 
erated RH cycling tests may be an effective way of simulating 
the automotive fuel cell duty cycles in order to observe the fail- 
ure mechanisms and to monitor the durability and the changes 
in the material properties of the membrane [7,21]. 


' Nafion is a registered trademark of E.I. DuPont De Nemours & Co. 
2 GORE-SELECT is a registered trademark of W.L. Gore & Associates Inc. 


In this work, we study the mechanical response of the mem- 
branes through numerical simulations. We employ a mechanics 
based model, that includes temperature and humidity dependent 
material properties, to determine the evolution of the stresses in 
the membrane due to the hydration—dehydration cycles at a fixed 
temperature (85°C) corresponding to the maximum operating 
temperature of the fuel cell. The influence of the cell assembly 
design in terms of the clamping method is investigated (Fig. 1). 
In addition, since the swelling of the membrane induces the 
compressive stresses and subsequent residual tension, swelling 
anisotropy of the membrane will be investigated. 

In the following, we review the theory for isotropic elasto- 
plasticity, followed by the definition of the geometry of the unit 
cell with boundary and loading conditions used in the numer- 
ical model. The material properties of the membrane used in 
this model are based on data obtained from our tensile tests 
of the membrane reported previously [17]. The details of the 
hygro-thermal loading cycles and the simulation of swelling 
strains during these cycles are studied. In Section 3, we intro- 
duce a series of assumed anisotropic swelling characteristics and 
evaluate their effects on the stresses. 


2. Theory for isotropic elasto-plasticity 


Here, we outline the approach used to incorporate hygro- 
thermal effects into the finite element simulation including 
isotropic plasticity. This is an extension of previous work, where 
first the linear-elastic [16], and later the linear-elastic, perfectly 
plastic behavior [14] of the membrane were considered. An 
uncoupled theory is assumed, for which the additional temper- 
ature changes brought about by the plastic strain are neglected. 

We assume that the total strain tensor, ¢;;, can be written as 
the sum: 


pl | oT S 
Eij = ef; + ej; + 8); + Ei (1) 


A ; l 
where Ejj are the elastic strain components, e}; are the plas- 


tic strain components, and ep and es are the temperature 
and swelling induced strains, respectively. Assuming isotropic 
behavior, the thermal strains resulting from a change in temper- 
ature of an unconstrained isotropic medium are given by 


1 ifi=j 


S _ z i ii = 
ay OE a eee 


(2) 


where a is the linear coefficient of thermal expansion cc, 
To and T are the reference and current temperature (°C), respec- 
tively, and 6; is the Kronecker delta. 

Similarly, the swelling strains caused by moisture uptake are 
given by 


e} = e} (T, H), (3) 


where T is the current temperature (°C) and H is the relative 
humidity (%). For an isotropic material, £;;(T, H) = e(T, H)ô;j. In 
the current work, ¢°(T, H) is defined as a temperature and humid- 
ity dependent polynomial function based on the experimental 
data [17]. 
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Fig. 1. The geometry of the unit cell used in the simulations. The mechanical boundary conditions are noted in the figure using symmetry assumptions. Two additional 
boundary conditions are imposed by subjecting the top of the stack to either constant displacement or constant force (not shown). 


Assuming a linear response within the elastic region, 
isotropic Hooke’s law is used to determine the stress tensor, 
Gij, and we have 


E 
a FA- 20) 


where Ekk = Exx + Eyy + Ez, E is Young’s modulus, and v is Pois- 
son’s ratio. Based on our previously published experimental data 
[17], Young’s modulus can be defined as a function of tem- 
perature and humidity, E = E(T, H), and we assume Poisson’s 
ratio remains constant. A generalized plane strain condition is 
imposed in the simulations: 


[vdizeq + (1 — 2v)ef, (4) 


Exy = Eyx = yz = zy =O and yy = Constant. (5) 


For the inelastic response, incompressible plastic deforma- 
tion is assumed, with the rate-independent plastic flow according 
to the von Mises yield function (J2-flow theory) [22]: 


3 
Foi =y 5 Su Si — 00, (6) 


where oj are the components of the (true) stress tensor, oo the 
yield strength, and Sj; are the components of the deviatoric stress 


tensor defined by 
1 
Sij = oij — 3 kK: (7) 
We further assume that the material exhibits isotropic hard- 
ening, thus the yield strength depends on the plastic strain as 
well as temperature and humidity, i.e.: 


oo = o0(é", T, H), (8) 


where ZP! is the current equivalent plastic strain: 


R 2 l 
gil = 5 dej dey. (9) 


In other words, this model considers isotropic expansion of 
the yield surface with plastic strain. According to the von Mises 
yield criterion, yield occurs when 


f(o;;, &') = 0, (10) 


and the material deforms elastically for f(ojj)<0. The Mises 
flow theory predicts that the plastic strain increment tensor is 
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proportional to the deviatoric stress tensor [22]: 


a1) 


a! = ggr! Fi, e) agr! Si 
K do, ij 2 00 
The equivalent plastic strain increment dé?! is determined 

from the condition that the yield criterion (Eq. (10)) must be 

satisfied during the strain hardening. 


3. Numerical model 
3.1. Assumptions 


A previously developed two-dimensional finite element 
model [14,16] is adapted for the current simulations. The model 
includes the following simplifying assumptions: 


(1) Simplified temperature profile with no heat generation. 

(2) The electrodes are integrated into the gas diffusion layer 
(GDL) to form a GDE instead of considering them as a 
separate layer. 

(3) For the GDE and the bipolar plate, the deformation is linear- 
elastic, with no swelling, while the membrane is allowed to 
deform plastically with isotropic hardening and swelling. 

(4) All material properties are isotropic including the thermal 
expansions. However, isotropic and anisotropic swelling of 
the membrane is investigated. 

(5) After the initial hygro-thermal loading to reach 85 °C, four 
humidity cycles are applied at this temperature, followed by 
hygro-thermal unloading. 

(6) During humidity cycling, a linear humidity gradient through 
the membrane is imposed from cathode to anode, such that 
cathode humidity is cycled from 30% RH to 95% RH while 
the anode is kept at 30% RH. 


3.2. Geometry and boundary conditions 


The gas flow channels in the graphite plates are used to trans- 
port hydrogen to the anode side of the membrane electrode 
assembly (MEA) and air to the cathode, while the product water 
vapor is carried away. In previous studies [14,16], we selected 
two graphite plate alignments to study the effects of the geome- 
try of the gas flow channels: (1) aligned and (2) alternating. We 
found that gas channel alignment primarily affects the location 
of the maximum stresses rather than the magnitude. Thus, for 
simplicity, in this study we will only investigate the aligned case 
in our simulations (see Fig. 1). 

We used the commercial software ABAQUS 6.4 [23] with 
four-node, generalized plane strain temperature—displacement 
coupled elements (CPEG4T) for the simulations. For compar- 
ison purposes, eight-node elements with reduced integration 
(CPEG8RT) were also tried. Since no significant difference was 
found in the results, in the following, we will present the results 
from the four-node elements. A total of 21,200 elements were 
used in the model and 9 x 80 elements are used for the thickness 
and length of the membrane, respectively. 


The thickness of the bipolar plates in the model is 11 mm, and 
the depth and the width of the gas channels in the plates are 1 mm 
each (Fig. 1). The membrane thickness is 50 um (corresponding 
to the thickness of Nafion® 112 membrane at ambient condi- 
tions). As mentioned previously, the catalyst layer is included in 
the gas diffusion layer, forming a gas diffusion electrode (GDE) 
with a thickness of 100 um. 

Two different clamping methods are investigated: (1) fixed 
force, corresponding to the case where the fuel cell stack is 
equipped with springs to control the clamping force, and (2) fixed 
displacement where the graphite plates are clamped preventing 
any net deformation in the thickness direction. In the former 
case, a constant pressure (1 MPa) is applied on the surface of 
upper graphite plate. In the latter case, a fixed displacement is 
calculated by applying a 1 MPa load to an otherwise unloaded 
model, and the resulting displacement at the upper boundary is 
fixed throughout the analysis. 

Referring to the unit cell in Fig. 1, symmetric boundary con- 
ditions, u,=0O at the bottom edge, and uy=0 on right edge 
of the unit cell are applied. In order to have uniform dis- 
placement at the left side of the membrane, a linear constraint 
Ajud! + Azu = 0 is applied on the left edge, where A; = 1 
and Az =—1, Nj is the node at the bottom left end and M2 are 
the remaining nodes on the left side. 


3.3. Material properties 


Previous work has shown that the mechanical response of 
the membrane is highly nonlinear [14]. It is therefore important 
to implement, as much as possible, the true material properties. 
The material properties for the graphite plates are set to that 
of commercial graphite and for the carbon paper the proper- 
ties are obtained from TORAY® TGP-H-030 [24]. It is assumed 
that these materials have linear thermo-elastic behavior and do 
not swell in the presence of moisture. The physical properties 
of the membrane and other components used in this study are 
shown in Table 1. Physical properties for the membrane were 
adopted from the DuPont™ Nafion® 112 (perfluorosulfonic 
acid) membranes data sheet [25] and listed in Table 1. For the 
mechanical properties of the membrane, we use the results of 
our experimental research, summarized in Tables 2 and 3 [17]. 
In the experimental study, tensile tests of the Nafion® 112 mem- 
brane films were conducted in an environmental chamber within 
a temperature range from 25 to 85°C, and a relative humid- 
ity range from 30% RH to 90% RH to obtain true stress-true 
strain curves at each temperature—humidity point. The details of 
the experiments are discussed elsewhere [17]. Based on these 
results, linear-elastic, plastic behavior with isotropic hardening 
and temperature and humidity dependent material properties 
is used for the membrane. Young’s modulus and the propor- 
tional limit (assumed to correspond to the yield strength) of the 
membrane are defined at four temperatures and four humidi- 
ties. Even though a slight anisotropy in the material properties 
was observed experimentally [17], in this study we assume, for 
simplicity, that the material properties are isotropic. 

The isotropic hardening used in the finite element pro- 
gram, ABAQUS [23], is defined by the initial yield strength, 
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Table 1 

The physical properties of the materials used in the finite element analysis [24,25] 

Units k(WmK~!) p (kg m~?) E (MPa) v a (1076 K7) Cp (J kg K7!) 
Bipolar plates (graphite) 95.0 1800 10,000 0.25 5.0 750 

GDE (carbon paper) 0.3 400 10,000 0.25 —0.8 500 
Membrane (Nafion® 112) 0.259 2000 Table 2 0.25 123.0 1050 


Table 2 
Young’s modulus at various temperatures and humidities for Nafion® 112 mem- 
brane [17] 


Relative humidity (%) 


30 50 70 90 
Young’s modulus? (MPa) 
T=25°C 197 192 132 121 
T=45°C 161 137 103 70 
T=65°C 148 117 92 63 
T=85°C 121 85 59 46 


à Linear interpolation between given data points is used during the loading 
and unloading. 


o(T, H; = 0), where the plastic strain is assumed to be zero, 
and two additional stress points at o!(T, H; P! =0.05), and o° (T, 
H; &P! =0.25). This is done for each temperature and humidity 
combination by using the available stress—strain curves of the 
material (Fig. 2). 

To incorporate the experimental data for the swelling strain 
(Fig. 3) [17] into the FE model, a third degree polynomial is used 
for the measured swelling strains at four humidities (30-90%) 
for each temperature. This way, the swelling strain can be written 


Table 3 

The data points of the true stress—plastic strain curve of Nafion® 112 mem- 
brane given at various temperatures and humidities [17] for (a) yield strength, 
o(eP! =0), (b) stresses at plastic strain 0.05, o! (eP! =0.05), and (c) stresses at 
plastic strain 0.25, o? (eP! =0.25) 


Relative humidity (%) 


30 50 70 90 
Yield strength (MPa) 
(a) P =0 
T=25°C 6.76 6.51 5.66 4.20 
T=45°C 5.67 5.21 5.01 3.32 
T=65°C 5.14 4.58 4.16 2.98 
T=85°C 3.61 3.44 3.08 2.20 
(b) e! =0.05 
T=25°C 7.16 6.61 6.22 5.11 
T=45°C 5.70 5.72 5.43 3.69 
T=65°C 5.30 4.77 4.36 3.33 
T=85°C 4.16 3.62 3.16 2.26 
(c) P! =0.25 
T=25°C 9.71 9.26 8.65 8.88 
T=45°C 7.31 7.34 7.48 6.18 
T=65°C 6.55 5.92 5.73 5.78 
T=85°C 5.04 4.28 4.22 4.31 


à Linear interpolation between given data points is used during the loading 
and unloading. 


Fig. 2. Schematic representation of the temperature and humidity dependent 
stress—strain curves with hardening. The values of stresses at each combination 
of temperature and humidity (7;, Hj) are given in Table 3(a-—c) for the plastic 
strain values, ¢?! =0, 0.05 and 0.25. 


as a polynomial function of humidity and temperature: 


4 
SC, H) = > Cyr iat, (12) 
i, j=l 


where Cj (i, j= 1. ..4) are the constants of the polynomial and 
given in Table 4. A reference state of 30% RH is chosen, where 
the swelling strain is assumed zero for all temperatures (Fig. 3). 
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Fig. 3. Experimental data (markers) for the swelling expansion in the membrane 
and the polynomial curve fit (solid lines) to these data points as a function of 
humidity and temperature, plotted for four constant temperatures for Nafion® 
112 [17]. 
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Table 4 
Constants for the swelling strain polynomial defined in Eq. (12) 
i=l i=2 i=3 i=4 

Cy 
j=l 2.994 x 10712 —5.221 x 10710 3.574 x 1078 —6.832 x 1077 
j=2 —4.303 x 107!0 7.361 x 1078 —5.166 x 1076 1.003 x 10-4 
j=3 2.163 x 1078 —3.566 x 1076 2.564 x 1074 —5.067 x 107° 
j=4 —5.402 x 1078 2.012 x 1075 —2.007 x 10-3 4.355 x 107? 


3.4. Simulation of hygro-thermal loading 


We simulate a fuel cell duty cycle by first defining the “initial 
conditions,” followed by conducting a loading sequence that can 
be summarized in four phases. 

The initial conditions (corresponding to the assembly con- 
ditions) for the zero stress-state are defined as all components 
of the cell stack being set to room temperature 20°C, and rela- 
tive humidity 30%. Since we do not have experimental data for 
temperatures below 25°C, the material properties below these 
values are assumed to be constant. The subsequent four load 
phases are: 


e Phase 1. The clamping condition (fixed force or fixed dis- 
placement as described above) is applied as described in 
Section 3.2. This corresponds to assembling the stack. 

e Phase 2. Hygro-thermal loading is applied by linearly increas- 
ing the relative humidity and temperature of the assembly 
simultaneously. A simplified temperature profile is imposed 
according to the assumption that the cathode/GDE and 
anode/GDE interface temperatures are at 85 and 86°C, 
respectively (Fig. 4). The bipolar plate is assumed to have a 
constant temperature at its midplane, 80°C, which is the cell 
operating temperature. These values correspond to a fuel cell 
operating at 80°C under typical conditions based on results 
from [2]. The humidity is modeled either as a constant RH 
value through the membrane “uniform loading” or as a gra- 
dient through the thickness “gradient loading.” In order to 


Graphite 


Cathode 
Side 
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30 95 80 90 
Relative Humidity (%) Temperature (°C) 


Fig. 4. The relative humidity and the temperature profiles in the membrane 
during the hydration—dehydration cycles at maximum thermal load. For the 
temperature profile of the unit cell, a constant anode/GDE temperature of 85 °C 
and cathode/GDE temperature of 86 °C is assumed [14,16], and a linear humidity 
gradient from cathode (fully humidified) through anode is imposed. 


simulate the humidity gradient distribution in the membrane, 
a linear variation is imposed in the thickness direction of the 
membrane such that the anode remains at the initial humidity 
value (30% RH), while the maximum humidity (95% RH) is 
reached at the cathode side (Fig. 4). These parameters roughly 
model the real-life operating conditions of a humidified mem- 
brane [2]. We note that the PEM is assumed to be the only 
component to swell in the presence of humidity. 

e Phase 3. Hydration—dehydration cycle is achieved by linearly 
decreasing and increasing the relative humidity between the 
initial (30% RH) and the hydrated (95% RH) conditions at 
the cathode, while the temperature is fixed at the maximum 
values (85/86 °C). This cycle is repeated four times (Fig. 5) 
(we will later see that a “cyclic steady state” is achieved, and 
further cycling will not provide additional information, within 
the context of this model). 

e Phase 4. Unloading is achieved by linearly decreasing the 
relative humidity and temperature simultaneously, back to the 
initial values, but leaving the clamping conditions imposed. 


The described sequence aims to simulate a typical accelerated 
test conducted to screen the performance and the durability of 
fuel cells [7,15]. This hygro-thermal loading cycle is a simplifi- 
cation of real-life operating condition that corresponds to taking 
the cell from start-up to a hydrated state and then back to the 
initial state. It is assumed that the cell remains at high tem- 
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Fig. 5. The relative humidity and temperature cycle in the membrane used for 
the analysis. Since there is a gradient for both temperature and humidity as given 
in Fig. 4, given here are the maximum values the membrane reaches: maximum 
humidity (95% RH) at cathode side and maximum temperature (86 °C) at anode 
side. 
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perature during the rapid hydration—dehydration cycles, since 
cooling occurs at a much slower rate than drying in an operating 
fuel cell. 


3.5. Simulating eigenstrain: swelling and thermal 
expansion 


The strains in the membrane due to the humidity and tem- 
perature changes are calculated from the assumed thermal and 
swelling strains (Eqs. (2) and (3)), where the coefficient of ther- 
mal expansion is given in Table 1. The user subroutine UEXPAN 
[23] is used to define the strain increments during the hygro- 
thermal loading/unloading in the simulations. 


3.5.1. Volumetric swelling strain 
The volumetric swelling strain, ev, is defined as 


ey = > = (1+ ŠA ted + es) 1, (13) 
where ¢;; are the components of the swelling strains in the normal 
directions (i=x, y or z) and dV and Vo are the volume change 
and the initial volume, respectively. 

Although swelling strains in the transverse and the machine 
direction exhibits a slight anisotropy [17], the swelling expan- 
sion behavior of Nafion® membrane is assumed isotropic 
[25]. At 85°C and 95% RH the swelling strain in Nafion® 
112 membrane is approximately 12% for both in-plane 
directions. Therefore, assuming isotropic material we have 
Exx = Eyy = Ex, = 0.12, where €,,, €yy are the in-plane (1, 2) and £37 
is the out-of-plane (3) swelling strains. Finally, the volumetric 
strain according to Eq. (13) is 
Ey = = = (1.12) — 1 = 0.40. (14) 

Vo 

In our simulations, we will first consider this assumed 
isotropic swelling followed by a study of anisotropic swelling 
effects on the mechanical response. 


3.5.2. Modeling of isotropic swelling 

When isotropic swelling is considered, the swelling strain 
is defined by a polynomial as a function of humidity and tem- 
perature, e5(7, H) according to Eq. (12). The swelling strain 
increment, AeS(T, H), due to a temperature—humidity incre- 
ment, is given by 
Ta, Hn) — &so(Tn—1, Hn-1), (15) 


1S0 


(A$ =$ 


1sO7n 1S0 


where n is the current loading increment and T, (°C) and H, (% 
RH) are the corresponding temperature and relative humidity 
values, respectively, at that increment. The subscript n — 1 indi- 
cates the condition at the previous increment and the subscript 
“iso” indicates isotropic swelling strain. The total swelling strain 
during the hygro-thermal loading/unloading is then 


N 
ego, H) = X Ae.) (16) 


n=1 
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Fig. 6. The evolution of the swelling strains during the first hygro-thermal cycle, 
from (30% RH at 20 °C) to the hydrated state (95% RH at 85 °C), and during the 
humidity unloading from hydrated state to the dry state (30% RH at 85°C) at 
constant temperature. The curves are plotted by using the swelling strain poly- 
nomial, e$ (T, H), whose constants are obtained by an interpolation among the 
experimental data (markers) obtained at given a set of temperature and humidity. 


where Nis the total number of loading increments. The evolution 
of the swelling strain during the loading and unloading is shown 
in Fig. 6. Due to the nonlinear dependence on temperature and 
humidity, separate paths are taken during loading and unloading, 
but the initial and end points are same (by definition). 


3.5.3. Modeling of anisotropic swelling 
We make two assumptions for simulating anisotropic 
swelling expansion of the membrane: 


(1) The total volume change of the fully hydrated anisotropic 
membrane due to the swelling is identical to the isotropic 
case (i.e. the maximum a = a = 0.40). This corre- 
sponds to requiring the same total water uptake. 

(2) The swelling expansions in the in-plane directions (x, y) are 
always equal (e$, = eS, Æ 2. ). 


Therefore, the anisotropy is introduced by increasing the 
swelling strain in the thickness direction (z), and calculating the 
corresponding swelling strains in the in-plane (x, y) directions 
using the constant volumetric swelling strain assumption. Three 
levels of anisotropy are investigated, ranging from the isotropic 
case to the “fully anisotropic” case (where the membrane swells 
only in the thickness direction). The swelling strain components 
calculated for each case are shown in Table 5. 

To investigate the influence of swelling anisotropy, we mod- 
ify the swelling strain increments for the isotropic case, by 
multiplying with the anisotropy ratio defined as 


; ; , : S, 
Anisotropic swelling strain — Eaniso 


Anisotropy ratio = 


Isotropic swelling strain a 
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Table 5 
The values of the total swelling strains in the membrane in each direction 
for the isotropic and anisotropic cases used in the analysis (anisotropy ratios, 


er fee (i =x, y) are given for each case, which are used to modify the 
subroutine) 
Isotropy Anisotropy 
Swelling strains (mm mm!) at maximum humidity (95% RH) load 
Directions 12-12 15-10 27-5 40-0 
In-plane: x 0.12 0.10 0.05 0 
In-plane: y (=x) 0.12 0.10 0.05 0 
Out-of-plane: z 0.12 0.157 0.27 0.40 
ee fe. 1 1.308 2.250 3.333 
7S 
oo al Ona 1 0.833 0.416 0 
where ee represents the anisotropic swelling strain in i- 


direction. For example, the swelling strain increments for the 
“15-10” case (swelling strain of 0.157 in the thickness direction 
and 0.10 in the in-plane directions) is 


S,z 
z Enis 0.157 
Ae! W= = A€iso = O12 dei = 1.308(Aéiso),(18) 


for thickness direction (z), and 


A 15-10 _ oe AEn = 0.10 -o = 0.833(Ag; ) (19) 
f= Ss Eiso = 0.12 Eiso = U. (Aéiso), 


1SO 


for in-plane directions (x, y), where (Ae'°) is given in Eq. (12). 
Other cases follow in a similar manner. The anisotropy ratios for 
the cases considered are given in Table 5. The thermal strains 
are assumed to be isotropic and the same for all cases. 

The resulting thermal and swelling strains — calculated 
through the subroutine UEXPAN - at the cathode are depicted 
in Fig. 7, in both in-plane and thickness directions for the first 
cycle. The “12-12” case represents the isotropic swelling prop- 
erties (Nafion® 112 membrane). The other curves represent the 


Isotropic 
(12-12) 


0.4 : 
Out-of-Plane 
sees In-Plane 

Ẹ Anisotropic 
£ (40-0) 
E 
E 
£ 
§ 0.2} ] 
” Anisotropic ees 
= (15-10) (27-5) 
3 
3 
”n 


0 D 
30%RH — 20 °C 


95%RH — 85 °C 30%RH — 85 °C 
Relative Humidity - Temperature during Cycle 


Fig. 7. The evolution of swelling strains during the first hygro-thermal cycle for 
the isotropic and three anisotropic cases, in the out-of-plane (thickness) (solid 
lines), and in-plane directions (dashed lines). The thermal strain is the same for 
all cases and plotted for comparison. 


three levels of anisotropy considered. Since the cathode and 
anode sides of the membrane are subjected to different rela- 
tive humidity values due to the assumed humidity gradient, the 
largest swelling strain occurs at the cathode side of the mem- 
brane, and there is no swelling at the anode side. We note that 
these strains drive the system. 


4. Results and discussion 
4.1. Isotropic swelling 


In a previous study [14] we found that, for the geometry con- 
sidered, it is sufficient to focus on two locations in the membrane 
to study the maximum and minimum range and magnitude of 
the stresses: left end and right end, see Fig. 1. Therefore, the two 
end-points will primarily be used to investigate the evolution of 
the stresses and plastic deformation in the following. We note 
that the left end is the midpoint of the gas channel groove and 
the right end is the midpoint of the land (Fig. 1). 

First, we investigate how the cyclic humidity loading affects 
the mechanical response of the membrane for two cases: 
(1) “uniform loading,” using a uniform humidity distribution 
through the membrane during cycling, and (2) “gradient load- 
ing,” using a humidity gradient from cathode to anode during 
cycling (Fig. 4). Uniform loading results in uniform in-plane 
stresses, whereas the gradient loading results in a stress gra- 
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Fig. 8. Effect of the humidity loading conditions in the membrane on the evo- 
lution of in-plane stress at the “left end” (A) and the “right end” (B) of the 
membrane (Fig. 1). The results shown are for two cases: uniform humidity and 
linear humidity gradient in the membrane. In the latter case maximum humidity 
is imposed at the cathode side of the membrane. 
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Fig. 9. Distribution of in-plane stresses in the membrane at maximum (95% RH at 85°C) and minimum (30% RH at 85°C) hygro-thermal load, and the residual 


stresses after unloading (30% RH at 20°C) (for the fixed displacement case). 


dient. Consequently, for the latter case, the in-plane stress is 
illustrated by showing the stresses at both the cathode and the 
anode side (Fig. 8). The gradient loading leads to similar stress 
levels at the cathode side of the membrane as those of the uni- 
formly loaded case but the stress at the anode side remains low, 
due to the low humidity there. This results in a significant stress 
gradient associated with gradient loading. Such a stress gra- 
dient can become associated with bending of the membranes 
if delamination occurs between the PEM and the GDE, possi- 
bly accelerating the failure evolution of the membrane. In the 
remainder of this study, the gradient humidity loading will be 
considered. 

The in-plane stress distribution in the membrane at maximum 
and minimum humidity loads and after unloading is depicted in 
Fig. 9 for gradient loading. The effect of the humidity gradient 
can be seen explicitly from these plots. The stresses at the cath- 
ode side of the membrane are larger than those of the anode side. 
Large compressive stresses are induced at maximum humidity 
load, and stresses are in tension at the minimum humidity load. 
After unloading, when both humidity and temperature decreases 
back to initial state, residual stresses up to 10 MPa are developed 
through the left of the membrane. Thus, the extreme maximum 
and minimum values of the stresses are localized at the left and 
right end of the membrane. 

We showed in Fig. 7 that the contribution of the thermal strain 
is very small compared to the swelling strains. Thus, the in-plane 
stresses are mostly driven by the swelling strain. Furthermore, 
the inelastic strain (in this study the plastic (permanent) strain) is 
an important contributor to fatigue induced failures in cyclically 


loaded structures [18]. It is therefore of interest to study the 
individual strain components during the cyclic loading. 

The contribution of each strain component in Eq. (1) (elastic, 
plastic, thermal and swelling) toward the total strain is plotted for 
the first cycle at the left and right ends of the membrane for the 
fixed displacement case (Fig. 10).? The initial hygro-thermal 
expansion is constrained by the surrounding structures (bipo- 
lar plates and GDEs), resulting in compressive elastic strains 
at both ends of the unit cell. The plastic strains remain zero 
until the yield limit is reached. When yielding starts, compres- 
sive in-plane plastic strains are induced due to the constraints at 
the ends. In the thickness direction, the resulting plastic strains 
are positive (due to incompressible deformation). Upon unload- 
ing, the elastic strains decrease, and go through a sign reversal 
eventually leading to reverse yielding. 

The evolution of the magnitude of the equivalent plastic strain 
= ij ij 
ing (added plastic strain) and reverse yielding (subtracted plastic 
strain) are observed for each cycle (Fig. 11C). We also note that 
plastic strains are different at different locations in the membrane 
due to the different constraints (Figs. 10 and 11C). 

The in-plane stress, oyy, and the out-of-plane stress, ozz, 
at the cathode side (Fig. 11A and B) both alternate between 
their respective upper and lower bounds, thus reaching a “cyclic 
steady state.” These stress extrema are determined by the loading 


4/(2/ 3)ePl Pl at the cathode side shows that forward yield- 


3 The elastic strains are not necessarily zero at the start of this cycle due to the 
initial clamping conditions imposed. 
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Fig. 10. (A and B) Contribution of thermal, swelling, plastic and elastic strains 
to the total strains in two directions; z (out-of-plane) and x (in-plane) for one 
humidity cycle, for the fixed displacement case. Due to the generalized plane 
strain assumption and the imposed boundary conditions, the strains in y-direction 
is equal to those of in x-direction. 


(temperature and humidity) and the yield strength of the mem- 
brane. Three important trends can be observed in Fig. 11A and B: 
(i) the magnitude of the out-of-plane stress is always smaller than 
the in-plane stress, (ii) at the left end, the in-plane stresses cycle 
between tensile and compressive loads, whereas the out-of-plane 
stresses cycle between two compressive values, and (iii) at the 
right end, all stress components cycle between two compressive 
values. Repeated fluctuation between tensile and compressive 
stresses, accompanied by cyclic yielding, is a more severe load 
case for fatigue compared to cyclic compressive stresses, and 
may be the cause of mechanical failure. We note that at the right 
end (midpoint of the land), the in-plane stresses remain com- 
pressive during cycling. Alternating between two compressive 
stress levels tends to be less likely to cause mechanical fatigue 
failure. 


4.2. Anisotropic swelling 


In this section, the effects of anisotropic swelling strain 
are investigated. The three anisotropic cases considered are 
designated; “10-15”, “27-5” and “40-0” (fully anisotropic) 
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Fig. 11. Evolution of (A) in-plane, (B) out-of-plane stresses, and (C) the magni- 
tude of equivalent plastic strain during the hygro-thermal cycles at three points of 
membrane’s cathode side for the fixed displacement case and humidity gradient 
loading. 


(see Section 3.5.3 for definitions and Table 5 for values). With 
the in-plane stress exhibiting the largest (cyclic) tensile stresses, 
as discussed above, we will limit the following discussion to the 
in-plane stress component, for brevity. 

The evolution of the in-plane stress at the two ends of the 
membrane is shown in Fig. 12 and corresponding plastic strain 
in Fig. 13 for the cases of anisotropic swelling considered. As 
can be seen from Fig. 12, increasing anisotropy leads to lower 
stresses at both ends. For the fully anisotropic case (40-0), the left 
end remains elastic during the cycling (Fig. 13A). For this case, 
the major contribution to the stresses comes from the clamping 
forces. For all cases with anisotropy, the right end experiences 
compressive stress (Fig. 12B) and yields only during the initial 
and final step (Fig. 13B) as opposed to the isotropic case, which 
exhibits cyclic yielding. 

The magnitude of the equivalent plastic strain at the left end of 
the model (the midpoint of the groove) is larger than at the right 
end (the midpoint of the land), and decreases with increasing 
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Fig. 12. The evolution of in-plane stresses for the fixed displacement case at (A) 
the left end and (B) the right end of the cathode side of the membrane plotted 
for different anisotropies. 
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Fig. 13. The evolution of equivalent plastic strain magnitude for the fixed dis- 
placement case at (A) the left end and (B) the right end of the cathode side of the 
membrane, plotted for different anisotropies (constant plastic strain corresponds 
to elastic response). 
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Fig. 14. The in-plane stresses, obtained after hydration and dehydration, respec- 
tively, during one hydration cycle plotted as a function of distance along the 
membrane for isotropic and anisotropic cases for fixed displacement. The arrows 
indicate the stress amplitude, the difference between the maximum and minimum 
stress, the membrane reaches during the cycle. 


anisotropy (Fig. 13). The membrane swells more in the in-plane 
direction with decreasing anisotropy and therefore the equivalent 
plastic strain magnitude at the left end increases. Constraints 
at the right end of the membrane prevent the expansion in all 
directions leading to a hydrostatic state of stress, which inhibits 
yielding according to the Mises yield criteria (Eqs. (6)—(10)). 

For the fully anisotropic case (40-0), the lack of swelling in 
the in-plane direction reduces the overall magnitude of stress 
resulting in vanishing plastic strain at the left end. No plastic 
strain develops at the left end. However, at the right end, large 
swelling of the membrane in the thickness direction leads to 
plastic deformation in the absence of hydrostatic stress. 

The in-plane stress distributions along the membrane at max- 
imum humidity load (95% RH) and after unloading (30% RH) 
are plotted for one cycle (the fourth cycle) in Fig. 14. In addition, 
the residual in-plane stresses at the cathode after cycling (end of 
the cyclic loading, 30% RH—20 °C) are plotted (for all cases) 
as a function of distance along the membrane, in Fig. 15. From 
Figs. 14 and 15 we see that the tensile stresses are highest at 
the left end of the model. This is due to the large plastic defor- 
mation occurring at these locations (Fig. 13A). At the right end 
(midpoint of the land), the smaller amount of yielding results 
in smaller residual stresses. For the fully anisotropic case, as 
discussed previously, there is no residual in-plane tensile stress 
for the left side of the model and most of the stresses come from 
the clamping load. However, for the right side, constraints pre- 
vent the thickening of the membrane, causing large compressive 
in-plane stresses. 


4.3. Effect of membrane thickness 
In order to examine the influence of the membrane thickness 


on the in-plane stresses, hypothetical membranes of thickness 
20 pm (with the same material properties as used above), are 
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Fig. 15. The distribution of residual in-plane stress along the membrane obtained 
at the end of the five cycles (after hygro-thermal unloading from hydrated state 
to the initial state; 30% RH—20 °C) for the isotropic and anisotropic cases, and 
fixed displacement. 


studied for the cases of isotropic (12-12) and fully anisotropic 
(40-0) swelling. In the previous subsections we have shown that 
the left and right ends of the model (midpoints of the groove and 
land) have the highest magnitudes of stresses during cycling 
(also see [14]). Thus, these stress amplitudes will be used to 
investigate the effect of membrane thickness. The in-plane stress 
amplitudes and the maximum stresses at the left and right ends 
of this hypothetical membrane (thickness: 20 um) are given in 
Fig. 16A and B for the two cases studied, compared with the 
original membrane of thickness 50 wm. 

For isotropic swelling, in-plane stress amplitudes, Ao = 


Bee oma % (Fig. 16A) and maximum in-plane stress, 
oKH—30% (Fig. 16B) both decrease for the thinner membrane, 


especially at the left end of the membrane. For fully anisotropic 
swelling, the thickness has little effect on the in-plane stresses. 
This stress reduction may seem counter intuitive (i.e. it seems 
a thinner membrane should give larger stresses), however since 
this is a strain-controlled problem, the thickness of the mem- 
brane plays a secondary role when determining the stresses. 
In summary, even though the membrane thickness affects the 
stresses, the change in swelling of the membrane is a more 
important factor governing the stresses. 


4.4. Effect of clamping conditions 


So far, the effects of anisotropy and thickness on the stresses 
and plastic deformation in the membrane have been examined 
for the fixed displacement clamping condition. As in the previ- 
ous subsection, stress amplitudes will be used to investigate the 
effect of the clamping conditions (fixed displacement and fixed 
force). Accordingly, the values of the stress amplitude, Ao = 


RH=30% _ ,RH=95% i RH=30% ; 
Onar = Onin , and the maximum stress, Oax , in 


the in-plane direction, for the left and right ends of the mem- 
brane, for both clamping methods are shown in Fig. 17 during 
the fourth cycle, for the isotropic and the anisotropic cases. 
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Fig. 16. (A) The in-plane stress amplitude (Ao, = gRE=30% 


(B) the maximum in-plane stress (GaN? ), at the left and right end of the 
membrane during the fourth cycle, for two thicknesses: 50 um and 20 um, for 
isotropic and fully anisotropic cases, and for fixed displacement case. 


RH=95% 
= Op. °) and 


The stress amplitude decreases significantly with increasing 
anisotropy (Fig. 17A) since the membrane swells less in the 
in-plane direction. Moreover, the in-plane stresses at the mini- 
mum humidity load, 30% RH (Fig. 17B), tend to decrease with 
increasing anisotropy. All maximum stresses are tensile except 
for the fully anisotropic case and for the right end of the fixed 
displacement case, where the maximum in-plane stress is com- 
pressive due to the absence of reverse yielding (Figs. 12 and 13). 
Therefore, this may be interpreted as an advantage of the fully 
anisotropic swelling since compressive stresses reduce the risk 
of propagating any existing flaws in the membrane. 

In previous studies, we investigated a case where the gas 
channels are not aligned. In general, the maximum stresses are 
similar to those of the aligned case and not shown here for 
brevity. However, in the alternating case, the large compressive 
stresses are in the midpoint of the model, instead of the right end 
[14]. 

In conclusion, increasing anisotropy results in a significant 
decrease in the stress amplitudes in the in-plane direction. Thus, 
it appears that increasing swelling anisotropy may be associated 
with reduced fatigue loading. 
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Fig. 17. (A) The in-plane stress amplitude (Ao, = oRH=30% — g =, and 


(B) the maximum in-plane stress (GhH= 02), at the left and right end of the 


membrane during the fourth cycle, for isotropic and anisotropic cases, and for 
fixed displacement and fixed force cases. 


5. Concluding remarks 


Hydration—dehydration cycles at constant temperature 
(~85 °C) are simulated in this numerical investigation to study 
the mechanical response of fuel cell proton exchange mem- 
branes, utilizing a unit cell approach. Both the case of uniform 
humidity loading and gradient humidity loading are consid- 
ered. When a humidity gradient over the membrane is assumed, 
the cathode is subjected to cyclic humidification but the anode 
remains at ambient humidity during the cycling, with the inten- 
tion of simulating conditions in an accelerated relative humidity 
(RH) cycling test. The mechanical model of the membrane 
assumes linear-elasticity with isotropic-hardening plasticity, and 
temperature—humidity dependent material properties based on 
our experimental studies of Nafion® membranes. 

Uniform humidity loading conditions result in a uniform 
stress distribution over the thickness of the membrane. Large 
in-plane tensile stresses are observed in the middle of the groove 
(“left end” of the model). For the case of a humidity gradient 
through the membrane thickness, a stress gradient is generated, 
leading to large stress amplitudes (where the stress amplitude 
denotes the difference between the maximum and minimum 


stress during one full cycle). Thus, within the framework of 
this study, mechanical failure is most likely to occur at the cath- 
ode side of the membrane. However, experimental evidence is 
needed to confirm these results. 

The results show that the in-plane swelling of the membrane 
due to water absorption is the primary driving force of the system 
for these types of load conditions, mainly causing large in-plane 
stresses. Our results show that swelling has a more profound 
effect on the stresses than either the clamping conditions or the 
membrane thickness. 

The influence of swelling anisotropy is also studied. To this 
end, a constant volume change due to swelling was assumed, 
where the decreased swelling strains in the in-plane directions 
correspond to an increased swelling in the thickness direc- 
tion. As expected, the in-plane stress amplitudes decrease with 
decreasing in-plane swelling, referred to as increasing swelling 
anisotropy. Since the in-plane stress is the major stress compo- 
nent in the membrane, and may cause the propagation of the 
minor flaws during cyclic loading, the parameters that affect 
the in-plane stress may have a direct influence on the fatigue 
behavior of the membrane. Furthermore, for fully anisotropic 
swelling, the membrane does not swell in the in-plane directions 
and the stresses never become tensile but remain compressive, a 
favorable outcome for the durability of the membrane. These 
findings suggest that membranes with low in-plane swelling 
strains during hydration—dehydration cycles may exhibit bet- 
ter performance during cell operation, which corresponds to the 
experimental results presented by Kolde et al. [13]. 

In conclusion, by designing the swelling anisotropy of the 
membrane, reduced stresses and therefore better fatigue resis- 
tance for PEM may be obtained that will ultimately improve the 
durability of the fuel cell. 
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